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RESPONSE OF DARK MATTER HALOS TO CONDENSATION OF BARYONS: 
COSMOLOGICAL SIMULATIONS AND IMPROVED ADIABATIC CONTRACTION MODEL 
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ABSTRACT 

The cooling of gas in the centers of dark matter halos is expected to lead to a more concentrated dark 
matter distribution. The response of dark matter to the condensation of baryons is usually calculated 
using the model of adiabatic contraction, which assumes spherical symmetry and circular orbits. In 
contrast, halos in the hierarchical structure formation scenarios grow via multiple violent mergers and 
accretion along filaments, and particle orbits in the halos are highly eccentric. We study the effects of the 
cooling of gas in the inner regions of halos using high-resolution cosmological simulations which include 
gas dynamics, radiative cooling, and star formation. We find that the dissipation of gas indeed increases 
the density of dark matter and steepens its radial profile in the inner regions of halos compared to the case 
without cooling. For the first time, we test the adiabatic contraction model in cosmological simulations 
and find that the standard model systematically overpredicts the increase of dark matter density in the 
inner 5% of the virial radius. We show that the model can be improved by a simple modification of 
the assumed invariant from M{r)r to M{f)r, where r and r are the current and orbit-averaged particle 
positions. This modification approximately accounts for orbital eccentricities of particles and reproduces 
simulation profiles to within 10 — 20%. We present analytical fitting functions that accurately describe 
the transformation of the dark matter profile in the modified model and can be used for interpretation 
of observations. 

Subject headings: cosmology: theory — dark matter : halos: structure — galaxies: formation — 
methods: numerical simulations 



1. INTRODUCTION 

During the past decade dissipationless cosmological sim- 
ulations have shown that the density distribution within 
virialized halos of different masses can be described by 
an approximately universal profile (Dubinski & Carlberg 
1991; Navarro et al. 1997; Moore et al. 1998). Affhough 
non-baryonic dark matter exceeds baryonic matter by a 
factor of ildm/^b ~ 6 on the average, the gravitational 
field in the central regions of galaxies is dominated by 
stars. In the hierarchical galaxy formation model the stars 
are formed in the condensations of cooling baryons in the 
halo center. As the baryons condense in the center, they 
pull the dark matter particles inward thereby increasing 
their density in the central region. 

The response of dark matter to baryonic infall has tradi- 
tionally been calculated using the model of adiabatic con- 
traction. Eggen, Lynden-Bell, & Sandage (1962) were the 
first to use adiabatic invariants of particle orbits to esti- 
mate the effect of a changing potential in a contracting 
proto-galaxy. Zeldovich et al. (1980) used the adiabatic 
invariant approach to calculate the contraction of lepton 
halos in response to the cooling of baryons and set con- 
straints on the annihilation cross-section of leptons. They 
also presented the first analytical expressions for adiabatic 
contraction (AC) for purely radial and circular orbits, as 
well as the numerical tests of such model. Barnes & White 
(1984) showed that the reaction of the particle orbits in 
spheroidal component to the slow growth of the disk in 
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their numerical experiments of disk and bulge evolution is 
indeed adiabatic and can be described by a simple model 
that assumes circular particle orbits and angular momen- 
tum conservation. 

The present standard form of the AC model was intro- 
duced and tested numerically by Blumcnthal et al. (1986, 
see also Ryden & Gunn 1987). This model assumes spher- 
ical symmetry, homologous contraction'^, circular parti- 
cle orbits, and conservation of the angular momentum: 
M{r)r = const, where M{r) is the total mass enclosed 
within radius r. With these assumptions, the final dark 
matter distribution is calculated given the initial mass pro- 
files Mdm(?'), M\,{r) and final baryon profile Mb(r/): 

[Md^ir) + Mb(r)] r = [Afdm(r) + Mb(r/)] r/. (1) 

This model has been studied further by Ryden (1988, 1991) 
and Flores et al. (1993). It is routinely used in mass mod- 
eling of galaxies (e.g., Flores et al. 1993; Dalcanton et al. 
1997; Mo et al. 1998; Courteau & Rix 1999; van den Bosch 
2001; van den Bosch & Swaters 2001; Klypin et al. 2002; 
Seljak 2002) and clusters of galaxies (e.g., Treu & Koop- 
mans 2002). The effect of the contraction of the dark 
matter distribution is important for studying star forma- 
tion feedback on the centers of dark matter halos (Gnedin 
& Zhao 2002) and for comparing the abundance of dark 
matter halos and galaxies as a function of circular velocity 
(e.g., Gonzalez et al. 2000; Kochanek & White 2001). It is 
particularly important in calculations of the dark matter 
annihilation signal from the Galactic center (e.g., Gnedin 
& Primack 2004; Prada et al. 2004). 

Despite recent advances in numerical simulations, the 
model of adiabatic contraction has never been tested in a 
cosmological context. The tests performed to date (Jes- 
seit, Naab, & Burkert 2002) assume spherical symmetry 

* The halo can be imagined consisting of spherical shells which con- 
tract in radius but do not cross each other. 
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and consider only the growth of a central concentration in 
an isolated halo. The hierarchical formation of halos is, in 
general, considerably more complex than the simple pic- 
ture of quiescent cooling in a static spherical halo. Each 
halo is assembled via a series of mergers of smaller halos, 
with the cooling of gas and contraction of dark matter oc- 
curring separately in every progenitor. The gas can be 
re-heated by shocks during mergers and during accretion 
along the surrounding filaments. Also, some objects may 
undergo dissipationless mergers after the gas is exhausted 
or the cooling time becomes too long. It was argued that 
dissipationless evolution erases the effect of gas cooling on 
the DM distribution (Gao et al. 2004). 

In this paper we consider the effect of dissipation on the 
dark matter distribution in high-resolution cosmological 
simulations. We also present the first test of the AC model 
in the self-consistent simulations of hierarchical structure 
formation and propose a simple modification which de- 
scribes numerical results more accurately. 

2. COSMOLOGICAL SIMULATIONS 

We analyze high-resolution cosmological simulations of 
eight group and cluster-sized and one galaxy-sized systems 
in a flat ACDM model: rij„ = 1 - f^A = 0.3, 17b = 0.043, 
h = 0.7 and erg = 0.9. The simulations are performed with 
the Adaptive Refinement Tree (ART) 7V-body-|-gasdynamics 
code (Kravtsov 1999; Kravtsov, Klypin, & Hoffman 2002), 
an Eulerian code that uses adaptive refinement in space 
and time and (non-adaptive) refinement in mass to achieve 
the high dynamic range needed to resolve the halo struc- 
ture. 

The cluster simulations have a peak resolution of « 
2.44/1-1 ]^p(, g^j^j j3]y[ particle mass of 2.7 x IQ^h^^ M© 
with only a region of ^ lOh^^ Mpc around each clus- 
ter adaptively refined. We analyze each cluster at a late 
epoch (0 < z < 0.43), when it appears most relaxed. This 
minimizes the noise introduced by substructure on the 
azimuthally-averaged mass profiles. The virial masses^ of 
the clusters range from « lO^/i-^ M© to 3 x lO^^h''^ Mq. 

The galaxy formation simulation follows the early (z > 
4) evolution of a galaxy that becomes a Milky Way-sized 
object at z = in a periodic box of 6h^^ Mpc. The simu- 
lation is stopped at z ~ 3.3 due to limited computational 
resources. At z = 4, the galaxy already contains a large 
fraction of its final mass: « 2 x 10^^h~^ M© within 30h~^ 
kpc. The DM particle mass is 9.18 x lQ^h~^ M© and the 
peak resolution of the simulation is 183h~^ comoving pc. 
This simulation is presented in Kravtsov (2003), where 
more details can be found. 

For each halo, we analyze two sets of simulations which 
start from the same initial conditions but include differ- 
ent physical processes. The first set of simulations follows 
the dynamics of gas "adiabatically" , i.e. without radiative 
cooling. The second set of simulations (hereafter CSF) 
includes star formation, metal enrichment and thermal 
supernovae feedback, metallicity- and density-dependent 
cooling, and heating due to the extragalactic UV back- 
ground. Star formation in the cluster simulations is imple- 
mented using the standard Kennicutt's law and is allowed 

^ We define the virial radius, rvir, as the radius enclosing an average 
density of 180 times the mean density of the Universe at the analyzed 
epoch. 
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Fig. 1. — Mass profile of one of the clusters as a function of 
physical radius. The solid and dotted lines show the profiles of dark 
matter and baryons (stars+gas) in the adiabatic (thin) and cooling 
(thick) runs, respectively. The dashed curve shows the prediction of 
the standard adiabatic contraction model, while dot-dashed curve 
shows the improved model. The profiles are truncated at four res- 
olution elements of the simulation. Top panel: relative mass differ- 
ence between the adiabatic contraction model and the DM profile in 
the CSF simulation. The dashed line is prediction of the standard 
AC model, while dot-dashed line shows our modified model. 



to proceed in regions with temperature T < lO^'K and gas 
density n > 0.1 cm"^. In the galaxy formation run, the 
star formation rate is proportional to the gas density and 
stars are allowed to form at densities n > 50 cm"'^. The 
difference in star formation prescriptions in galaxy and 
cluster simulations, accounts for the difference in spatial 
resolution. The prescription used in the galaxy formation 
run is more appropriate when applied at the scale of tens 
of parsecs (see Kravtsov 2003, for discussion). 

To identify dark matter halos we use a variant of the 
Bound Density Maxima algorithm (Klypin et al. 1999). 
Dark matter particles in the high-resolution region of the 
simulation are assigned a local density calculated by using 
24-particle SPH kernel. We identify local density peaks on 
a scale of lOOh^^ kpc and analyze the density distribution 
and velocities of the surrounding particles to test whether 
a given peak corresponds to a gravitationally bound ob- 
ject. In this study we only consider host halos: those that 
do not lie within a larger virialized halo. We identify the 
center of each halo with the position of the DM particle 
with the highest local density. Based on the convergence 
studies for the ART code (Klypin et al. 2001; Tasitsiomi 
et al. 2004), we truncate the dark matter profiles at the 
inner radius 4Aa;inin, where Axmin is the smallest cell size: 
2.44/i~i and 0.183h~^ comoving kpc in the cluster and 
galaxy formation runs, respectively. 

3. EFFECTS OF COOLING ON MATTER DISTRIBUTION 
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Fig. 2. — Density profile of the cluster shown in Figure 1, with 
the same line types. In order to emphasize the differences at small 
radii, we plot the combination p{r)r^ which is roughly constant 
for isothermal distributions. Physical radius is shown in h~^ kpc 
(bottom axis) and as a fraction of the virial radius, r^i,- {top axis). 



We start by comparing the spherically averaged distri- 
bution of baryons and dark matter in the adiabatic and 
CSF simulations. The comparison reveals the effect of 
cooling and star formation because the two simulation for 
each object have the same initial conditions. Note that 
the cluster simulations likely suffer from the "overcooling" 
problem: the fraction of gas in the cold phase is about a 
factor of ~ 2 — 3 higher than suggested by observations 
(e.g., Balogh et al. 2001). The effect of cooling on mass 
distribution is thus likely overestimated. For the purposes 
of the present study, however, this is acceptable. In fact, 
the larger effect of cooling allows us to emphasize the dif- 
ference between the simulations and the model. 

The mass and density profiles of a representative cluster 
are shown in Figures 1 and 2. These figures show that 
cooling leads to an increase in the dark matter density 
within r < 50/i~^ kpc or r/rvir < 0.04 (see also Tissera 
& Dominguez-Tenreiro 1998). The mass profile is affected 
substantially at r < O.lrvir and the change increases with 
decreasing radius. At larger radii the average mass profile 
is not sensitive to baryon dissipation and the differences 
between the adiabatic and CSF runs are the result of the 
slightly different location of massive substructures. 

Figure 3 shows the density profiles in the galaxy for- 
mation run at 2; = 4. Qualitatively, the effect of cooling 
is similar to that seen in the cluster simulations. In this 
case, however, DM density in the CSF run is enhanced 
within a larger radius, r/rvir ^ 0.1. Also, baryons domi- 
nate the total density at r/r^iv ^ 0.03 in the galaxy simu- 
lation, while in the cluster simulation they dominate only 
at r/rvir ^ 0.01. The difference is due to the consider- 
ably higher fraction of cold (T < 10^ K) gas in the galaxy 
run: 80% at z = 4 versus ~ 0.2 — 0.3 in the cluster sim- 
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Fig. 3. — Density profile in the galaxy formation run at ^ = 
a function of physical radius. Lines types are as in Figure 1. 
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ulations. These fractions reflect a difference in densities, 
temperatures, and cooling times of the gas in clusters and 
high-redshift galaxies. Note that there are differences in 
the mix of gas and stars. At z < 0.5, most of the baryon 
mass within the central regions of clusters is in the stel- 
lar component of the cD galaxy. In the galaxy run, on 
the other hand, more than 80% percent of baryons in the 
dense central disk are still in the gaseous form. The cluster 
and galaxy formation simulations thus probe qualitatively 
different regimes of the evolution of central baryon con- 
densation. 

4. TESTING THE ADIABATIC CONTRACTION MODEL 

4.1. Standard Model 

In this section we test the standard prescription for adi- 
abatic contraction, given by equation (1). In order to 
calculate the model prediction for the final dark matter 
distribution in the CSF run, we use the mass profile of 
baryons Mb(ry) from this run and the mass profiles of 
dark matter and baryons from the adiabatic simulation at 
the same epoch. We consider adiabatic profiles as the ini- 
tial profiles for the model, Mdm(r) and Mb(r). We then 
use equation (1) to predict the DM distribution in the CSF 
run and compare the model prediction to the actual dark 
matter profile in simulation. 

The prediction of the standard AC model is shown in 
Figures 1, 2, and 3 by the dashed lines. The fractional 
deviations of the model prediction from the simulation for 
all analyzed systems are shown in Figure 4. The standard 
model predicts the overall mass enhancement but system- 
atically overestimates its magnitude in the inner regions, 
at r/rvir ^0.1. This effect has already been noticed in 
the original study by Blumenthal et al. (1986). Possible 
causes of the discrepancy could be (1) non-spherical mass 
distribution and substructure; (2) simultaneous evolution 
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of the dark matter and baryonic components; and (3) the 
assumption of circular orbits. In the next section wc inves- 
tigate whether the model can be improved by accounting 
for orbital eccentricities. 

4.2. Modified model 

The orbits of particles in dark matter halos in simula- 
tions are highly eccentric (e.g., Ghigna et al. 1998). The 
combination M{r)r for such orbits varies with the orbital 
phase and is not an adiabatic invariant. The conserved 
quantities for eccentric orbits are the angular momentum, 
J, and the radial action, 



1 



Vr dr, 



(2) 



where Vr is the radial velocity, and rp and Tq are the peri- 
and apo-center, respectively. For non-crossing spherical 
shells, the radial velocity can be expressed using the first 
integrals of motion C (e.g., Ryden & Gunn 1987; Gnedin 
& Ostriker 1999) as 

/2GM(r) J2 \i/' 

Vr = — - — + C . (3) 

\ r r^ / 

The parameters C change when shells cross, although the 
sum over all shells is constant. In the case of purely ra- 
dial orbits in a self-similar potential, we have J = and 
rp = 0, so that /^ oc M{ra)ra, and therefore the combi- 
nation M{ra)ra is conserved during a slow change of the 
potential (Blumenthal et al. 1986). The potentials of dark 
matter halos are, in general, not self-similar while the an- 
gular momentum is not zero, and the invariants cannot be 
reduced to a simple combination of r and M{r) that would 
be useful for predicting the final dark matter profile. 

It is interesting to check, though, if a combination that 
depends on the average radius along the orbit, f, rather 
than the instantaneous value r or the maximum r^, is con- 
served better than M{r)r or M{ra)ra. The orbit-averaged 
radius is 

- 2 p dr .^. 

^^TJT r—, (4) 

-^r Jrp Vr 

where Tr is the radial period. 

As a first check, we consider isochrone potential, <I>iso oc 
-[1 + {l + (r/rs)2)i/2]-i (Binney & Tremaine 1987), for 
which the radial action is known analytically as a function 
of E and J. We choose the angular momentum that corre- 
sponds to realistic eccentricities (e ~ 0.7): J — Jc{E)/^/3, 
where Jc{E) is the angular momentum of a circular orbit of 
energy E. We then check whether the combination M{f)f 
is proportional to I^ for all orbits. The ratio M{f)f/I^ is 
constant at small radii {r <^ rg) and large radii {r ^ rg) 
but varies in between by 50%. An alternative, the apoccn- 
ter ratio M{ra)ra/Ir exhibits considerably larger variation 
between the asymptotes, by about 300%. Thus the com- 
bination M{f)f is a better proxy to the invariant than 
M{ra)ra. 

In order to study more realistic orbits with the ener- 
gies and angular momenta relevant for cosmological simu- 
lations, we use a separate set of high-resolution collision- 
less simulations of three Milky Way-sized halos and one 
Virgo cluster-sized halo. Galaxy-sized halos have about 
one million DM particles within their virial radius (see 
Kravtsov et al. 2004, for details), while the cluster-sized 
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Fig. 4. — Fractional diflferences of the mass profiles predicted 
by the adiabatic contraction models and the simulation profiles for 
eight clusters (top eight panels) and one galaxy formation run (the 
bottom panel). Dashed lines correspond to the standard model and 
dot-dashed lines show our modified model, eq. (6). The cluster 
in the top panel is experiencing a merger event with a comparable 
mass cluster, which can be seen as an excursion of the profile at 
r ~ 0.2rvir. 



halo has approximately eight million particles (Tasitsiomi 
et al. 2004). All runs have spatial resolution < lO'^^rvir. 
We approximate the potential of each halo with a spher- 
ically symmetric NEW profile (Navarro et al. 1997) pa- 
rameterized by the mass and maximum circular velocity 
measured in the simulations. We then integrate the or- 
bits in this potential starting with the particle positions 
and velocities given by the simulations and calculate the 
radial action (eq. 2) numerically. 

The orbits have a distribution of eccentricities that is 
very similar to the isotropic distribution in analytical po- 
tentials studied by van den Bosch et al. (1999): the 20%, 
50%, and 80% quartiles are e = 0.39,0.61,0.79, respec- 
tively. Given such a wide eccentricity distribution, the 
orbit-averaged radius f varies for particles at a given cur- 
rent radius r depending on the orbital phase. Nevertheless, 
the mean relation at 10""^ < r/r^-„ < 1 can be described 
by a power law function 

X = Ax'" , x = r/rvii, (5) 

with small variations in the parameters A and w from halo 
to halo and from epoch to epoch. The mean values are 
A w 0.85 ± 0.05 and w « 0.8 ± 0.02, which we use as our 
fiducial parameters. This power-law dependence reflects 
typical energy and eccentricity distributions of particles in 
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cold dark matter halos. At a; < 0.44 the average radius 
is larger than the current radius, while at a; > 0.44 it is 
smaller. 

The distribution of eccentricities at a current radius r 
leads also to a distribution of each combination M{r)r, 
M{ra)ra, and M{f)f. None of these combinations is pro- 
portional to the invariant I^ for individual orbits. What 
we are looking for is a proxy to the invariant for an en- 
semble of orbits representing a spherical shell. Such a 
proxy would allow us to calculate the transformation of 
the dark matter profile imagined of consisting of spherical 
shells. Therefore, instead of individual particle orbits we 
look at the average ratios M{ra)ra/l'^ and M{f)f/I^ in 
radial bins. Across all bins the average ratio M{ra)ra/I^ 
varies by a factor 3.4, while the ratio M{f)f/I^ varies only 
by a factor 2.2. Therefore, our insight gained from analytic 
orbits in isochrone potential is still valid for isotropic or- 
bits in NFW potential: the combination M{f)f is a better 
proxy to the invariant. Moreover, we have found that the 
mixed combination M{f)r/I^ varies even somewhat less, 
by a factor of 2 across all radii. This latter combination 
thus is our best proxy for the radial action. 

Motivated by these considerations, we propose a modi- 
fied adiabatic contraction model based on conservation of 
the product of the current radius and the mass enclosed 
within the orbit-averaged radius: 



M{r)r = const. 



(6) 



Using equation (5) we compare our gasdynamic simula- 
tions with this modified model. 

Figures 1-4 show that the modified model provides a 
more accurate description of the simulation results than 
the standard model for most of the objects. Although 
there are still deviations of the model and simulation pro- 
files, there is no systematic over- or under-prediction of 
the mass for all objects. Typical deviations are < 10%. 
Note that for the high-redshift galaxy disk both the stan- 
dard and modified AC models overpredict the mass profile 
significantly, although the modified model is still closer to 
the simulated profile. The discrepancy is large in the in- 
ner kiloparsec where the most of the mass is in the thin 
gaseous disk. The adiabatic contraction model cannot be 
applied in such regime and further numerical simulations 
are needed to investigate the dark matter distribution in 
the galactic center. 

An interesting question is why the adiabatic contrac- 
tion approximation works as well as it does? Although the 
model assumes that gas cooling affects DM distribution in 
the final object adiabatically, the particles experience con- 
traction in separate unconnected halos and undergo some 
degree of violent relaxation in subsequent mergers. To 
check the evolution of individual particle orbits we have 
compared their physical distance r to the center of the 
main progenitor and quantities M{r)r, M(f)r at a number 
of epochs from z = to z = 4 for one of the clusters. There 
is a substantial scatter in the relation of these quantities 
between the present and higher redshifts. However, some 
interesting trends for the averages of M{r)r and M(f)r 
can be observed. First of all, we find that both radii and 
products of radii and mass evolve as the object grows hi- 
erarchically, especially during major mergers. However, 
between mergers, during the periods when a substantial 
amount of gas cools, M{r)r and M{f)r seem to be con- 
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Fig. 5. — The effect of gas cooling on the mass profile at different 
epochs for the cluster shown in Figures 1 and 2. Thin and thick lines 
show the profiles at 2 = 4 and z = 0.18, respectively: dotted lines 

— profiles for baryons (gas+stars) in the CSF run; dot-dashed lines 

— profiles of DM in the adiabatic run; solid lines — DM profiles in 
the CSF run. 



served. The latter combination is conserved better than 
the former, which may explain why our modification to 
the AC model is successful. 

We find also that most (~ 70%) particles within the 
central 25h~^ kpc of the cluster center at z = come from 
a single progenitor at z = 1. This progenitor is not the 
most massive but has the highest central density. When 
it merges with a more massive progenitor at z « 0.5 to 
form the final system, its particles are most tightly bound 
and end up dominating the mass in the central region of 
the merger remnant. This is consistent with other merger 
simulations which show that particles from the progenitor 
with the densest central region dominate the inner regions 
of the remnant (Boylan-Kolchin & Ma 2004). 

5. HOW STEEP IS THE CENTRAL DARK MATTER 
PROFILE? 

Dark matter distribution in the central regions of galax- 
ies and clusters is important for testing the CDM paradigm 
and interpreting the observations. It is generally believed 
that dissipation by baryons would steepen the dark mat- 
ter profile (Barnes & White 1984; Blumenthal et al. 1986; 
Ryden & Gunn 1987). However, Loeb & Peebles (2003) 
and Gao et al. (2004, hereafter G04) recently suggested 
that after an early epoch of cooling and rapid star for- 
mation, subsequent dissipationless mergers would erase 
the cooling-induced central concentration of dark matter. 
They put forth a hypothesis that the NFW profile is a 
dynamical attractor, in the sense that remnants of dissi- 
pationless mergers are driven to the central profile with 
the cusp p(r) oc r~^, even if their progenitors have steeper 
DM profiles. As a supporting argument G04 use obser- 
vation that the density in the inner regions of halos in 
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dissipationless simulations is approximately constant af- 
ter the initial period of rapid mass growth. The attractor 
hypothesis together with observation that stars dominate 
gravitationally in the centers of galaxies leads to the con- 
clusion that the dark matter profile in galaxies which ex- 
perienced dissipationless mergers should be shallower than 

The results of our simulations do not support this con- 
clusion. Given the importance of the problem, it is worth 
discussing the differences between the simulations and anal- 
ysis of G04. Figure 5 compares the mass profile of the 
cluster shown in Fig. 1 at z = with its most massive 
progenitor aX z — A. We first verify that in adiabatic sim- 
ulations the DM mass profile in the inner w lO/i^^ kpc is 
approximately constant from z — A to the present epoch, 
in agreement with analyses of Fukushige & Makino (2001) 
and G04. Within the physical radius r = 10/i~^ kpc the 
enclosed mass is M{r) « 5 x 10^^ h^^ Mq at both epochs. 

The evolution is very different in the run with cooling 
and star formation. By z — 4, a, considerable stellar mass 
(~ 3 X 10" h~'^ Mq) has formed within 10 h'^ kpc of the 
center of the main cluster progenitor. However, the baryon 
mass within central 10 h~^ kpc continues to grow and in- 
creases by a factor of ten between z = 4 and z = 0.2. 
Approximately 50% and 70% of those stars form at z < 1 
and z < 2, respectively. Note that the stars and cold gas 
in cluster cores are accumulated both due to direct cooling 
in the core and via accretion during mergers (Motl et al. 
2004). 

As can be seen in Figure 5, such substantial increase in 
the baryon mass leads to the increase of the DM density 
in the inner regions. The final dark matter mass within 
10 h^^ kpc is w 2 X 10^^ h~^ Mq, or a factor of four larger 
than the mass in the adiabatic simulation. Therefore, one 
of the major differences between our simulations and anal- 
ysis of G04 is that in the simulation the density of both 
baryons and dark matter increases at lower redshifts, while 
G04 assume that the density of DM decreases as the total 
mass distribution is driven to the NFW profile. Note also 
that G04 assume that star formation and cooling in the 
centers of massive ellipticals effectively stops at z ^ 2 — 3, 
while in our simulations cooling and star formation con- 
tinue at lower redshifts with more than half of the stars 
formed at z < 2. This can be a deficiency of simulations 
which, as we mentioned above, suffer from the overcool- 
ing problem. The fraction of baryons in the cold gas and 
stars within the virial radius of clusters at z = in our 
simulations is in the range ~ 0.3 — 0.4, at least a factor 
of two higher than observed for the systems of the mass 
range we consider (Lin, Mohr, & Stanford 2003). The 
overcooling and relatively late star formation is a generic 
problem of cosmological simulations and is hardly realis- 
tic. It likely indicates that some mechanism suppressing 
cooling is needed. 

If we follow G04 and assume that cooling ceases at late 
epochs, the question is then whether subsequent dissipa- 
tionless evolution erases the prior effect of cooling on the 
concentration of DM and drives the overall stellar+DM 
profile to the NFW form, as required by the attractor hy- 
pothesis. Several recent studies have considered the ef- 
fect of major mergers and dynamical effects of substruc- 
ture on the dark matter profiles in the inner regions of 



halos (Dekel et al. 2003a,b; El-Zant et al. 2003; Boylan- 
Kolchin & Ma 2004; Ma & Boylan-Kolchin 2004; Nipoti 
et al. 2004). Dissipation of the orbital energy of massive 
subhalos by dynamical friction can heat the dark matter 
particles of the host and reduce the DM density in the in- 
ner regions. On the other hand, DM particles of the sink- 
ing subhalos replace the host particles and the overall cen- 
tral density profile stays approximately constant or even 
becomes steeper, depending on internal structure, spatial 
distribution, and orbital parameters of subhalos (Ma & 
Boylan-Kolchin 2004). 

Major mergers can lead to a more drastic and violent 
rearrangement of matter in halos compared to the effects 
of substructure. Boylan-Kolchin & Ma (2004) present a 
set of merger simulations, which are most relevant for our 
discussion. They show that mergers of halos with constant 
density cores produce a remnant with a constant density 
core, but mergers of the cored and cuspy halos produce a 
cuspy remnant. Their analysis shows that the initial cusp 
is remarkably stable and the density distribution of the 
merger remnants retains memory of the density profiles 
of their progenitors. This is in good agreement with the 
analysis of merger experiments presented by Kazantzidis 
et al. (2004). The merger remnant of the DM halos with 
a steep inner density profile (p(r) oc r^^*) retains the 
initial slope of the inner cusp. The mergers of halos with 
embedded stellar disks also produce DM halos with the 
inner cusp not shallower than the initial. 

Finally, to test the effect of dissipationless mergers on 
the DM density profile and the validity of the AC model, 
wc repeated the cluster simulation shown in Figure 5 with 
the cooling turned off at z < 2. The evolution is thus dissi- 
pationless at low redshift. The cluster undergoes a number 
of minor mergers, as well as one approximately equal-mass 
merger at z ~ 0.5. The amount of cold gas and stars within 
the virial radius of this cluster at z = is « 14%, a factor 
of two smaller than in the original simulation and similar 
to the stellar fractions observed in galaxy clusters (e.g., 
Lin et al. 2003). We find that although the effect of con- 
traction in this case is smaller, as less baryons condense 
in the center, the dark matter profile is still steeper than 
in the adiabatic run and is well described by the adiabatic 
contraction model. Dissipationless mergers that the main 
cluster progenitor has undergone at z < 2 apparently have 
not erased the steepening of the profile due to the cooling. 

The attractor hypothesis is thus not supported both by 
dissipationless merger simulations and our simulation in 
which cooling is stopped at high redshift. The effects of 
dissipation on the DM distribution in the progenitors is 
retained, at least to a certain extent, in the density distri- 
bution of their descendant. The overall effect of merging 
is just to mix dark matter particles within the same dis- 
tribution, while baryon dissipation leads to a significant 
increase of the dark matter density. 

Finally we note that for the halos that harbor a su- 
permassive black hole the density distribution within its 
sphere of infiuence (typically 10 — 100 pc) is determined 
by the interaction between the black hole, stars, and DM 
(e.g., Gnedin & Primack 2004). 

6. DISCUSSION AND CONCLUSIONS 
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We have analyzed results of self-consistent cosmological 
simulations of eight galaxy clusters and one galactic halo 
with and without cooling and star formation. The compar- 
ison of adiabatic and CSF simulations shows that cooling 
increases the total density and the density of dark matter 
at r < O.lrvir- This agrees qualitatively with results of 
other recent simulations (Tissera & Dominguez-Tenreiro 
1998; Lewis et al. 2000; Pearce et al. 2000; Valdarnini 
2002). Note that at r > O.Olrvir modern dissipationless 
simulations have reached a robust convergence (e.g., Die- 
mand et al. 2004). We conclude therefore that further 
progress in making predictions for the DM distribution on 
small scales requires studying gas dissipation. The effect 
of cooling also needs to be taken into account when com- 
paring the simulations with observations. 

We have presented the first tests of the adiabatic con- 
traction model in self-consistent high-resolution cosmolog- 
ical simulations.^ We find that the standard AC model 
systematically overprcdicts the increase of the dark mat- 
ter density in the inner < 0.05rvir. We have shown that 
the model can be improved by a simple modification of 
the assumed conserved invariant from M{r)r to M{f)r, 
where r and f are the current and orbit-averaged particle 
positions. This modification approximately accounts for 
the eccentricity of particle orbits. Our improved model 
describes profiles in simulations considerably better than 
the standard model, with the average accuracy of 10—20%. 

Jesseit et al. (2002) have used controlled simulations of 
isolated spherical halos with a growing central concentra- 
tion of baryons to test the validity of the standard AC 
model. They find that generally the standard model de- 
scribes their simulation more accurately than we find in 
our tests against cosmological simulations. Possible causes 
of the discrepancy with our results are the differences of 
(i) orbital distributions for the isolated halos versus cos- 
mological halos formed by mergers; and (ii) formation his- 
tories and degree of violent relaxation. In particular, Jes- 
seit et al. (2002) find that the standard model significantly 
overpredicts the contraction effect when the growth of the 
central baryon concentration is rapid compared to the dy- 
namical time of the halo. 

In order to understand this discrepancy, we have run 
isolated simulations of a live NFW halo with an analyti- 
cal contracting disk using the publicly available tree code 
GADGET (Springel et al. 2001) with 10^ particles. In 
agreement with Jesseit et al. (2002), we find that both 
the standard and the modified models work well (within 
20%) in the regions where the baryon density is compara- 
ble to the final dark matter density. However, the modified 
model becomes progressively more accurate as the baryons 
dominate over the dark matter by a factor of 2-3. In this 
case the ratio of the baryon density to the initial DM den- 
sity is a factor of ~ 20 — 50, while the ratio of the final 
to initial DM density (compression ratio) is ~ 10. For 
example, the modified model should apply in the central 
parts of the Galaxy. From these experiments we conclude 

® After this paper was completed we have learned of a different study 
of the AC model in cosmological simulations (Gottbrath & Stein- 
metz, 2000 unpublished). These authors found that the AC model 
works adequately at the resolved scales. However, the resolution of 
the simulations used in their study is considerably lower than res- 
olution of our simulations, and the scales where we find significant 
discrepancy between AC model and simulations were not resolved. 



that the accuracy of our model is determined both by the 
orbital distributions and by the amount of compression. 

Given that dark matter halos assemble via mergers and 
violent relaxation, it is somewhat of a puzzle that the adi- 
abatic contraction model reproduces the results of simu- 
lations so well. The success of the model seems to im- 
ply that the effect of the central baryon condensation on 
the dark matter distribution is independent of the way in 
which this condensation is assembled. At the same time, 
it may simply be due to the fact that the central region is 
dominated by particles from a single densest progenitor. 
If the progenitor halo contracts in response to the cooling 
of baryons early on and then approximately preserves the 
shape of its inner density profile during subsequent merg- 
ers, as suggested by merger simulations (see § 5), the AC 
model applied to the final mass distribution is expected to 
work. 

In Appendix we provide analytical fitting functions that 
describe the contraction of an initial NFW profile in our 
modified model. These functions can aid in interpretation 
of observations of galaxy halos and clusters of galaxies. We 
show that the inner slope of the dark matter density profile 
7 is determined by the shape of the baryon profile (see 
eq. [A12]). For the specific cases of the exponential disk 
and Hernquist model, both the baryon and the contracted 
DM profiles have the asymptotic slope 7 = 1. 

Our results have several implications for the efforts to 
test predictions of the CDM model observationally. The 
test that received much attention in the last decade is 
the density distribution in the inner regions of galaxies 
and clusters. Extensive convergence studies have shown 
that modern highest-resolution dissipationless simulations 
agree in their predictions: the average logarithmic slope of 
the density profile at r = O.Olrvir is 7 « 1.3 with a substan- 
tial scatter of ±0.3 from object to object (Fukushige et al. 
2004; Tasitsiomi et al. 2004; Navarro et al. 2004; Reed et al. 
2004; Diemand et al. 2004). At the same time, despite a 
significant decrease in the smallest reliably resolved scale, 
the logarithmic slope continues to get shallower with de- 
creasing radius without reaching an asymptotic value. For 
the purposes of the present discussion, it suffices that dis- 
sipationless simulations have converged at the scales where 
the effects of dissipation become important. 

Observational measurements of the dark matter density 
distribution are notoriously difficult. Several approaches 
have been used but in each case opposite conclusions are 
reached by different researchers, often after analyzing the 
same data. The rotation curves of dark matter dominated 
dwarf and low-surface brightness galaxies tend to favor 
density profiles shallower than predicted by CDM (e.g., Si- 
mon et al. 2003; de Blok et al. 2003). The results, however, 
are sensitive to the resolution of rotation curves, presence 
of bulges and non-circular motions, and for many galaxies 
the profiles can be reconciled with theoretical expectations 
(e.g., Swaters et al. 2003; Rhee et al. 2004) (but see de 
Blok 2003). The analysis of density distribution for bright 
galaxies is complicated by the uncertain contribution of 
stars to the total mass profile (e.g., Treu & Koopmans 
2002; Mamon & Lokas 2004). Some analyses tend to fa- 
vor inner slopes shallower than predicted by CDM (e.g.. 
Gentile et al. 2004), but others deduce slopes of the inner 
profiles (at least marginally) consistent with predictions 
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(Treu & Koopmans 2002, 2004; Koopmans & Trcu 2003; 
Jimenez et al. 2003). 

The density distribution in galaxy clusters can, in prin- 
ciple, provide a cleaner test of the models because effects 
of "gastrophysics" on the DM distribution are expected to 
be smaller and simpler. The advances in lensing analyses. 
X-ray observations, and high-resolution spectroscopy al- 
low observers to obtain constraints on the DM distribution 
using a variety of techniques. Recent observational anal- 
yses using strong lensing (e.g., Tyson et al. 1998), high- 
resolution X-ray imaging (Ettori et al. 2002; Katayama & 
Hayashida 2004, but see David et al. 2001 and Arabad- 
jis et al. 2002), and lensing+velocity dispersion measure- 
ments for the central galaxy (Sand et al. 2002, 2004) seem 
to favor shallow (7 < 0.5 — 1) inner density distributions. 
If these observational results stand, it would be a major 
problem for the CDM because dissipation generally makes 
the discrepancy worse. 

Many of the systematic effects and validity of the key 
assumptions in such measurements are yet to be explored. 
For example, deviations from spherical symmetry in the 
mass distribution allows for the inner slopes of 7 > 1 
(Dalai & Keeton 2004; Bartelmann & Meneghetti 2004). 
Some of our results also call into question the assumptions 
used to derive observational constraints. For example, it 
is often assumed that at small scales the DM profile can 
be well approximated by a power law, p{r) ex r~'^ . It is 
also assumed that dissipation would steepen the profile 
predicted from dissipationless simulations while retaining 
its power law form. As we show in Appendix, however, the 
steepening of the profile due to cooling is in general scale- 
dependent. For realistic cases, the profiles of dark matter 
and baryons at r < O.Olrvir should be quite similar. These 
scales are exactly where the profiles of massive ellipticals 
and central cluster galaxies are probed by spectroscopic 
measurements in observations. If incorrect assumptions 
are made about the DM distribution, there is a danger of 
oversubtracting the contribution of baryons to the total 
profile. Stars dominate the inner region and only a small 
overestimate of the stellar mass-to-light ratio could lead 
to a much lower residual density of dark matter. 

The fact that the standard AC model overpredicts the 
effect of cooling on mass distribution means that obser- 
vational analyses that use the standard model (e.g., Treu 
& Koopmans 2002, 2004) provide somewhat less stringent 
limits on 7 than claimed. Our results for the disk galaxy at 
high redshift show that even the improved model can still 
overestimate the effect in the inner region of the gaseous 
disk. Extrapolation of the model to very small radii, as 
is often done in predictions of the dark matter annihi- 
lation signal, can therefore be dangerous. Further high- 
resolution gasdynamics simulations are needed to probe 
the effect of cooling in the centers of dark matter halos. 



was supported by the National Science Foundation under 
grants AST-0206216 and AST-0239759 to the University 
of Chicago, by NASA through grant NAG5-13274, and 
in part by grant by the Kavli Institute for Cosmological 
Physics (KICP), an NSF Physics Frontier Center, through 
grant NSF PHY-0114422. O.Y.G. is supported by the 
STScI Fellowship. D.N. is supported by the NASA Grad- 
uate Student Researchers Program and by NASA LTSA 
grant NAG5-7986. The simulations presented here were 
performed on the IBM SP4 system (copper) of the Na- 
tional Computational Science Alliance. 



We would like to thank Frank van den Bosch, Leon 
Koopmans, Joel Primack, David Weinberg, and Simon 
White for discussions and useful comments on the manuscript, 
and Andrew Zentner and Stelios Kazantzidis for discus- 
sions and analysis of the density profiles in the controlled 
merger experiments. We are grateful to Nicole Papa for 
careful reading of the manuscript. We would also like to 
acknowledge Bocage of San Saba Vineyards. This work 



ADIABATIC CONTRACTION OF DARK MATTER HALOS 



APPENDIX 



ANALYTICAL FITS FOR THE CONTRACTED DARK MATTER MASS PROFILE 

In many applications it can be useful to have an analytical formulae to estimate the compression of dark matter due to 
baryonic condensation. We find that a simple procedure described below provides an accurate fit for the compression of 
an initial NFW profile. 

We assume that the initial distributions of dark matter and baryons are both given by the NFW profile, Mi{r), with 
a concentration parameter c. Subsequently the baryons cool and form stars and their final profile is given by Mb(r/), 
with the ratio of the baryon-to-total mass at the virial radius Mb(rvir)/Mvir = /b- This ratio does not need to equal the 
universal baryon fraction and may deviate from it depending on the details of hierarchical formation and heating by the 
extragalactic UV flux. We consider two representative examples of the final baryon distribution: an exponential disk, 
which is appropriate for spiral galaxies, and a Hernquist model, which adequately describes the stellar profiles in elliptical 
galaxies. Both of these profiles are characterized by a scalelength r\^. 

With the usual assumption of spherical homologous contraction, our goal is to calculate the final radius rf of the dark 
matter particles initially located at r > rf. It is useful to define the contraction factor y = rj/r. The equation we need 
to solve for y is 

rM,{f) = [(1 - /b)M,(f) + Mb(f/)] rf. (Al) 

Let us divide both sides of the equation by Mi{r)rf and use dimensionless radius x = r/r^i-c and mass m = M /M^-^: 

1 1 . , ™b(Sf) 

- = J- - Jb H TTT-- 

y m^{x) 

By definition, there is no contraction at the virial radius: y(l) = 1. 

An important simplifying feature of the profiles we consider here (NFW, exponential disk, Hernquist model) is that at 
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X <$^ 1 the enclosed mass of all three profiles grows with radius as m{x) ex x 
at a; <IC 1 are: 
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where gc = [log(l + c) — c/(l + c)]" 



Since our fit for the orbit-averaged radius x is also a power-law {x — Ax^ , w — 0.8) 
the mass ratio in eq. (A2) is constant at a; ^ rb, c^^, and therefore y{x) is approximately constant. 

The finite value yo = y(0) is pivotal to constructing the fitting function for y{x). We obtain the following equation for 
j/o by retaining only the first terms of the Taylor expansion: 

l = l-/b + ay2-, (A4) 

where a = 2/b(l + ^b)^/(^b.9cC^) for the Hernquist model and a = .fh / (r'^gcC^) for the exponential disk. There is no 
analytical solution for this algebraic equation, but we can obtain a very accurate approximate solution as follows. We 
first consider the standard model w — 1, for which eq. (A4) becomes a cubic equation. The solution is 

Q'^' + ^X''-(q'^'-IX'\ Q^C-^X + tJ^- (A5) 
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There is only one real solution because Q is always positive. We then note that eq. (A4) with w = 0.8 differs from the 
cubic equation only at a > 1. We can therefore obtain a solution for yo by matching the solution of the cubic equation 
for small a with the correct asymptote for large a. We find that the following modification produces the right asymptotic 
solution for a > 1: 

l/(l + 2u.) / -. ^ l/(l + 2«,) . / , „ ^ l + 2tu , 

)i/2^^^ _fni/2_^\ n=L(}—A\ ^J_ (A6) 



Vi'L. 



1 

2a 



1 
2a 



l + 2wj (2a)2' 

Finally, we link the two asymptotes with a weight function that we empirically found to minimize the error of the fit: 

yo - y(i) e-2" + y(„) (1 - e-2a) . (A7) 

The relative error of this fit, compared with the direct numerical solution, is only 0.2% on the average for all a and /b. 
The asymptotes a ^ 1 and a ^ 1 are reproduced almost exactly, while the maximum error of 0.6% occurs at a ^ 1. 

Now we need to connect the boundary values 2/(0) = j/o a-nd y{l) = 1 by an interpolating function. This can be done 
very accurately with the following trick. Let us define an auxiliary function 
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Fig. A6. — Increase of the dark matter mass for the initial NFW models assuming the exponential disk for the final baryon profile. Solid 
lines are the numerical solution, while dashed lines are our fitting formulae. Five sets of lines show that our analytical fit (eq. [Alf]) reproduces 
the direct numerical solution for various combinations of the parameters c, /t,, and r^,, with a typical accuracy of better than a few percent. 



It satisfies the boundary conditions t(0, j/o) = 2/0 f^nd t(l, 1) = 1. Since y{x) varies slowly near the two bounds, we can use 
the function t(x, y) to calculate the asymptotic solutions: y{x ^ 1) « t{x, yo) and y{x « 1) ss t{x, 1). The two asymptotes 
can be linked by a smooth weight function: 

y(x) = t(x, yo) e-"- + t(x, 1) (l - e-"-) . (A9) 

The exponent b that minimizes the error of the fit can be found by an approximate Taylor expansion of eq. (A2) near 
a; = 0. We find 
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where n = 1 for the Hernquist model and n = 3 for the exponential disk. We have tested the accuracy of eq. (A9) against 
direct numerical solution for y[x) varying the parameters in the range 4 < c < 20, 0.02 < /b < 0.24, 0.01 < rb < 0.07, 
which should cover most of the cases of interest. We find an average relative error of 2% for the Hernquist model and 1% 
for the exponential disk. The maximum error is 7% and 4%, respectively, and occurs at a: ~ rb- 

We can now easily find the final dark matter profile at the contracted radius Xf = xy{x): Mdm.iixy) = Afdm,i(a^)- If 
desired, the profile can be re-mapped to the grid of initial radii x by interpolation. Let us define the compression function 
Fm as the ratio of the final dark matter mass to the initial mass at the same radius xy{x): 



Fnixy) 



Mdms{xy) Mdm,i(a;) 



(All) 



MdiaAxy) Mdm.i{xy) 

We see that the increase of the dark matter mass due to baryonic infall can be calculated from the initial dark matter 
profile M(3in,i(a;) given the transformation y{x). Since we can use an analytic NFW profile, the accuracy of calculating 
Fm{x) is similar to that for y{x). Figure A6 illustrates that eq. (All) provides an accurate fitting function at all radii 
for all values of three independent parameters (c, /b, ^b). 

To summarize, in order to calculate the DM mass profile after contraction analytically one needs to: 

(1) calculate the maximum compression value j/o using eqs. (A5-A7); 

(2) calculate the exponent h of the weight function using eq. (AlO); 

(3) calculate the interpolating function y{x) using eqs. (A8~A9); 

(4) calculate the increase of the enclosed dark matter mass using eq. (All). 

As mentioned above, the mass profiles of the NFW, exponential disk, and Hernquist models are similar at a; <C rb, c~^: 
m(x) (X x"^ . This coincidence leads to a finite enhancement factor for the dark matter mass and density in the central 
region: Afdm,f(0)/Afdma(0) = Pdm,f(0)/pdm,i(0) = y^^- Therefore, the central contraction factor y^ (eq. [A7]) provides 
a useful measure of the maximum enhancement due to baryonic condensation. The constant contraction factor means 
that the inner slope of the dark matter distribution after contraction is the same as before contraction. The radius at 
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FlG. A7. — The slope of the dark matter profile in cluster MS 2137-23. Dotted lines are the initial NFW profile with c = 5 and c = 3. Dashed 
and solid lines are the post-contraction profile for the final baryon distributions given by the Hernquist model and Jaffe model, respectively, 
with /b = 1.3 X 10~^ and r^ = 0.016. Dot-dashed lino shows the post-contraction profile in the case of the Jaffe model with c = 3, /b = 10~^ 
and rb = 0.016. 



which this slope is reached, however, depends on the halo concentration and the baryon scale length (see Figure A6). At 

intermediate radii the post-contraction slope is steeper. In general case, the slope of the final DM profile is different from 

the initial profile and depends on the shape of the baryon density profile. 

Since the value of the inner slope of the dark matter density profile is important for testing CDM models against 

observations, it is interesting to consider the change of the slope due to the condensation of baryons into a configuration 

with an arbitrary final density profile. Let the baryon density at a; <C 1 be Phix) oc x"" and therefore the mass 

m\,{x) ex x^~" . If i^ > 1, the contraction would not be self-similar with a; j ex a; in the center. Instead, the post-contraction 

radius would scale as some power of the initial radius: x^ ex x". Substituting this into equation (A2) and letting a; ^ 0, 

we find a = (1 -I- 2i(;)/(l + (3 — v)w). The post-contraction dark mass at the center would scale as ra(x) ex x^/", and 

therefore the density is 

1 -|- Iwv 
p{x)^x-\ j=-——. A12 

1 + 2w 

For v — I, we recover our previous result, 7 = 1, which means that after the contraction of NFW profile by baryons which 
form an exponential or Hernquist profile, the asymptotic slope at small scales remains the same. For larger i/, the slope 
7 becomes steeper than the initial slope. Equation (A12) shows (recall that w ~ 1) that for t^ ^ 1 — 2 the inner slope of 
the dark matter profile will be quite close to the slope of the baryonic profile. 

As an illustration, we consider the profiles with the parameters similar to those of the cluster MS 2137-23 described by 
Sand et al. (2002, 2004): the virial mass, radius, and concentration are Mvir ~ 10"'^^ M© and rvir ~ 2 Mpc, and c « 5, 
respectively. At the center of the cluster lies a cD galaxy, which can be well fit by the Jaffe model: pb oc a:~^(a; -I- rb)~^ 
with Tb w 33 kpc. Assuming the best fitting mass-to-light ratio M/Ly — 2.6 (Sand et al. 2004), the mass of the central 
galaxy is 1.2 x 10^^ Mq. In our notation, the relevant parameters are /b = 1.2 x 10"'^, r^, — 0.016. 

Figure A7 shows the slope of the post-contraction dark matter density profile appropriate for MS 2137-23 in the range 
10^'^ < X < 10~^ which is used in the Sand et al. analysis. If the baryon distribution was described by a Hernquist model 
{v = 1), the inner slope would tend to 7 = 1, as we discussed above. For the Jaffe model (v = 2), on the other hand, the 
slope steepens at a; < 10~^. According to equation (A12), the asymptotic inner slope is 7 w 1.6. This slope is reached 
faster if we consider a more massive galaxy (/b = 10^^) and a less concentrated cluster {c ~ 3). 
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